function [wave_coher,wave_phase,time_out,scale_out,coi_out,cross_wavelet]=wave_coherency(wave1,time1,scale1,wave2,time2,scale2,coi1,dt,dj);

%******************************************************************* WAVELET
%+
%   WAVE_COHERENCY
%
% PURPOSE:   Compute the wavelet coherency between two time series.
%
%
%
% INPUTS:
%
%    WAVE1
%
%
% OUTPUTS:
%
%    WAVE is the WAVELET transform of Y. This is a complex array
%    of dimensions (N,J+1). FLOAT(WAVE) gives the WAVELET amplitude,
%    ATAN(IMAGINARY(WAVE),FLOAT(WAVE)) gives the WAVELET phase.
%    The WAVELET power spectrum is ABS(WAVE)^2.
%
%
% OPTIONAL KEYWORD INPUTS:
%
%    DT = amount of time between each Y value, i.e. the sampling time.
%         If not input, then calculated from TIME1(2)-TIME1(1)
%
%    DJ = the spacing between discrete scales.
%         If not input, then calculated from SCALE1
%
%
% OPTIONAL KEYWORD OUTPUTS:
%
%    SCALE_OUT = the vector of scale indices, given by S0*2^(j*DJ), j=1...J
%            where J is the total # of scales.
%
%    COI_OUT = the vector of the cone-of-influence
%
%----------------------------------------------------------------------------
% Copyright (C) 1998-2000, Christopher Torrence
% This software may be used, copied, or redistributed as long as it is not
% sold and this copyright notice is reproduced on each copy made.  This
% routine is provided as is without any express or implied warranties whatsoever.
%
% Notice: Please acknowledge the use of the above software in any publications:
%    ``Wavelet software was provided by C. Torrence and G. Compo,
%      and is available at URL: http://paos.colorado.edu/research/wavelets/''.
%
% Reference: Torrence, C. and G. P. Compo, 1998: A Practical Guide to
%            Wavelet Analysis. <I>Bull. Amer. Meteor. Soc.</I>, 79, 61-78.
%
% Please send a copy of such publications to either C. Torrence or G. Compo:
%  Dr. Christopher Torrence               Dr. Gilbert P. Compo
%  Advanced Study Program                 NOAA/CIRES Climate Diagnostics Center
%  National Center for Atmos. Research    Campus Box 449
%  P.O. Box 3000                          University of Colorado
%  Boulder CO 80307--3000, USA.           Boulder CO 80309-0449, USA.
%  E-mail: torrence@ucar.edu              E-mail: gpc@cdc.noaa.gov
%----------------------------------------------------------------------------
% modif Charria Guillaume Dec 2006 (IDL -> MATLAB)
% Update: time and scale smoothing - Elodie Gutknecht
%****************************************************************** WAVELET



%*** find overlapping times
    if (min(time1) > min(time2))
    	time_start = min(time1);
    else
        time_start = min(time2);
    end
    if (max(time1) > max(time2))
        time_end = max(time2);
    else
        time_end = max(time1);
    end
        time1_start = find( time1 == min(time1(find(time1 >= time_start))) );
        time1_end = find( time1 == max(time1(find(time1 <= time_end))) );
	time2_start = find( time2 == min(time2(find(time2 >= time_start))) );
	time2_end = find( time2 == max(time2(find(time2 <= time_end))) );

     
%*** find overlapping scales
    if (min(scale1) > min(scale2))
    	scale_start = min(scale1);
    else
        scale_start = min(scale2);
    end
    if (max(scale1) > max(scale2))
        scale_end = max(scale2);
    else
        scale_end = max(scale1);
    end
        scale1_start = find( scale1 == min(scale1(find(scale1 >= scale_start))) );
	scale1_end = find( scale1 == max(scale1(find(scale1 <= scale_end))) );
	scale2_start = find( scale2 == min(scale2(find(scale2 >= scale_start))) );
	scale2_end = find( scale2 == max(scale2(find(scale2 <= scale_end))) );
    
%*** cross wavelet & individual wavelet power
	cross_wavelet = wave1(time1_start:time1_end,scale1_start:scale1_end).*conj(wave2(time2_start:time2_end,scale2_start:scale2_end));
	power1 = abs(wave1(time1_start:time1_end,scale1_start:scale1_end)).^2;
	power2 = abs(wave2(time2_start:time2_end,scale2_start:scale2_end)).^2;

    



switch nargin
case 7
    dt = time1(2) - time1(1);
    dj = log(scale1(2)/scale1(1))/log(2);
case 9
otherwise
    disp('Wrong number of input arguments')
    wave_coher=[];
    wave_phase=[];
    time_out=[];
    scale_out=[];
    cross_wavelet=[];
    coi_out=[];
    return
end


	
	ntime = time1_end - time1_start + 1;
	nj = scale1_end - scale1_start + 1;
	
	scale = scale1(scale1_start:scale1_end);

% % ---- Modified version by Elodie ----
%    	for j=1:nj; %*** time-smoothing
%     		nt = fix(fix(4*scale(j)/dt)/2*4) + 1;
%     		time_wavelet = ([1:nt] - nt/2).*dt./scale(j);
%     		wave_function = exp(-time_wavelet.^2./2.);   %*** Morlet
%     		wave_function = wave_function./sum(wave_function); % normalize
%     		nz = fix(nt/2);
%     		%zeros_var(1:nz) = complex(0,0);            
%     		%cross_wave_slice = [zeros_var,cross_wavelet(:,j)',zeros_var];
%                 %cross_wave_slice = conv(cross_wave_slice,wave_function); 
%                 cross_wave_slice = conv(cross_wavelet(:,j)',wave_function);
%       		size(cross_wave_slice);
%                 size(cross_wavelet);            
%                 nz+1;
%                 ntime+nz;
%                 cross_wavelet(:,j) = cross_wave_slice(nz+1:ntime+nz)';                       	
%                 %zeros_var = double(zeros_var);
%     		%power_slice = [zeros_var,power1(:,j)',zeros_var];
%     		%power_slice = conv(power_slice,wave_function);
%                 power_slice = conv(power1(:,j)',wave_function);
%     		power1(:,j) = power_slice(nz+1:ntime + nz)';
%                 %power_slice = [zeros_var,power2(:,j)',zeros_var];
%     		%power_slice = conv(power_slice,wave_function);
%                 power_slice = conv(power2(:,j)',wave_function);
%     		power2(:,j) = power_slice(nz+1:ntime + nz)';           
%     end %*** time-smoothing

%    % ---- Original version ----
%     	for j=1:nj  %*** time-smoothing
%     		nt = fix(4*scale(j)/dt)/2*4 + 1;
%     		time_wavelet = ([1:nt] - nt/2).*dt./scale(j);
%     		wave_function = exp(-time_wavelet.^2./2.);   %*** Morlet
%     		wave_function = wave_function./sum(wave_function); % normalize
%     		nz = nt/2;
%     		zeros_var(1:nz) = complex(0,0);
%            
%     		cross_wave_slice = [zeros_var,cross_wavelet(:,j)',zeros_var];
%                cross_wave_slice = conv(cross_wave_slice,wave_function); 
%       		cross_wavelet(:,j) = cross_wave_slice(nz+1:ntime+nz);
%     		zeros_var = double(zeros_var);
%     		power_slice = [zeros_var,power1(:,j)',zeros_var];
%     		power_slice = conv(power_slice,wave_function);
%     		power1(:,j) = power_slice(nz+1:ntime + nz);
%     		power_slice = [zeros_var,power2(:,j)',zeros_var];
%     		power_slice = conv(power_slice,wave_function);
%     		power2(:,j) = power_slice(nz+1:ntime + nz);
%         end %*** time-smoothing

  
    %*** normalize by scale
    scales(1:ntime,1:nj) = NaN;
    for iscale=1:ntime
        scales(iscale,:)=scale;
    end
    
  
	cross_wavelet = cross_wavelet./scales;
	power1 = power1./scales;
	power2 = power2./scales;

    nweights = fix(0.6/dj/2 + 0.5)*2 - 1;   % closest (smaller) odd integer
    weights(1:nweights)=1.;
    weights = weights./sum(weights);  % normalize
    size(weights);
   
% % ---- Original version ----    
%    for i=1:ntime % *** scale-smoothing
%            cross_wavelet(i,:) = conv(cross_wavelet(i,:),weights);
% 	    power1(i,:) = conv(power1(i,:),weights);
%            power2(i,:) = conv(power2(i,:),weights);
%    end %*** scale-smoothing
%   
% % ---- Modified version by Elodie ----
%     for i=1:ntime % *** scale-smoothing
%         nz = fix(nweights/2);
%         cross_wave_slice2 = conv(cross_wavelet(i,:)',weights);
%         cross_wavelet(i,:) = cross_wave_slice2(nz+1:nj+nz)';
%         cross_power_slice2 = conv(power1(i,:),weights);       
%          power1(i,:) = cross_power_slice2(nz+1:nj+nz);
%         cross_power2_slice2 = conv(power2(i,:),weights);
%         power2(i,:) = cross_power2_slice2(nz+1:nj+nz);
%     end %*** scale-smoothing

 
    wave_coher = (abs(cross_wavelet).^2)./(max(power1.*power2,1e-9));
	wave_phase = 180./pi.*atan2(imag(cross_wavelet),real(cross_wavelet));%degre
    %wave_phase = angle(cross_wavelet); % radian
    
	time_out = time1(time1_start:time1_end);
	scale_out = scale1(scale1_start:scale1_end);
    if (length(coi1)==length(time1))
        coi_out=coi1(time1_start:time1_end);
    end
    



